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Executive  Summary 


Probabilistic  engine  health  management  (PHM)  is  expected  to  be  a  go-forward  approach  for  the 
USAF  and  other  DoD  agencies  to  enable  dramatic  improvements  in  the  assessment  and 
management  of  military  assets.  As  a  result,  accurate  and  information-rich  probabilistic  lifing 
methods  are  essential  to  assess  the  benefits  of  technology  insertion  programs  for  PHM. 

The  objective  of  this  research  was  threefold:  1)  to  enhance  probabilistic  lifing  methods  to 
provide  new  information  regarding  the  sensitivity  of  probabilistic  lifing  estimates  with  respect  to 
non-destructive  inspection  methods,  material  parameters,  loading,  stress  levels,  etc.,  1)  to  address 
computational  efficiency  issues  such  that  application  of  the  enhancements  are  practical,  and  3)  to 
integrate  the  new  methods  with  existing  approaches  as  permissible  to  foster  adoption. 

Three  major  tasks  were  undertaken:  1)  develop  sensitivity  methods  for  lifing  scenarios  with 
simulated  inspections  processes,  2)  explore  new  complex  variable  methods  for  sensitivity 
analysis,  and  3)  enhance  the  score  function  method  to  address  new  scenarios. 

In  each  case  significant  research  was  accomplished  that  lead  to  several  enhancements.  In 
particular:  1)  a  new  probabilistic  sensitivity  method  was  developed  that  will  compute  the  partial 
derivative  of  the  probability-of-failure  with  respect  to  the  parameters  of  a  Probability-of- 
Detection  (POD)  curve  for  negligible  cost,  2)  new  complex-variable  sensitivity  methods  were 
developed  and  exploited  to  determine  shape  sensitivities  for  finite  element  models  (particularly 
with  respect  to  crack  size)  in  a  highly  accurate  manner  and  to  solve  initial  value  differential 
equations  with  extreme  accuracy,  and  3)  a  new  method  to  compute  probabilistic  sensitivities  with 
respect  to  bounds  of  truncated  distributions  was  developed  and  tested. 

Ancillary  objectives  were  also  accomplished.  In  particular,  four  students  were  supported  wholly 
or  in  part  in  obtaining  a  Masters  of  Science  in  Mechanical  Engineering,  4  conference  papers 
were  submitted  and  presented,  1  journal  paper  was  published,  2  journal  papers  are  under  review, 
and  4  Matlab-based  software  programs  were  provided  to  AFRL  for  technology  transition.  The 
software  programs  consisted  of  a)  a  generic  probabilistic  sensitivity  module,  b)  a  weight 
function-based  fatigue  lifing  code,  c)  an  inspection  sensitivity  program,  and  d)  complex-variable 
finite  element  code  for  shape  sensitivity  analysis. 
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1  Sensitivity  methods  probability-of-failure  estimates  with  respect  to 
POD  curve  parameters 

1.1  Introduction 

Nondestructive  evaluation  (NDE),  also  known  as  nondestructive  testing  (NDT)  or  nondestructive 
inspection  (NDI),  plays  a  vital  role  in  fracture  control  plans.  Methods  such  as  visual,  dye 
penetrant,  ultrasonics,  radiography,  and  eddy  current  are  among  the  common  inspection 
techniques  used  to  ensure  structural  integrity  [1],  The  type  of  inspection  and  the  times  of 
inspection  must  be  carefully  selected  to  ensure  safety  with  reasonable  cost.  There  are  a  number 
of  industries  that  have  a  long  history  of  application  of  NDE  methods  for  structural  integrity.  For 
example,  applications  include  nuclear  [2-4],  petroleum  [5],  aircraft  structures  [6-8],  gas  turbines 
[9-10],  and  offshore  structures  [11]  to  name  a  few.  Rummel  at  al.  [12]  provide  a  summary  of  a 
number  of  issues  related  to  the  application  of  NDE  methods  to  systems. 

The  inspection  process  is  simulated  using  a  Probability  of  Detection  (POD)  curve  that  quantifies 
the  probability  of  detecting  a  crack  as  a  function  of  the  crack  size.  A  probabilistic  lifing  analysis 
integrates  probability  distributions  of  crack  sizes,  load,  geometries,  and  material  properties  with 
the  inspection  processes  to  determine  the  probability-of-failure  with  and  without  inspection.  The 
purpose  is  to  ensure  a  sufficient  level  of  reliability  in  the  structure  and  to  optimize  the  inspection 
process  for  cost  and  maintenance-related  issues. 

However,  there  is  always  some  doubt  as  to  the  correct  values  for  the  parameters  of  a  POD  curve 
for  a  particular  field  application  due  to  the  mismatch  of  the  environment  for  development  of  the 
POD  curve  versus  its  field  application.  In  addition,  there  may  be  value  in  “what- if”  scenarios  to 
quickly  assess  the  value  of  changes  in  the  inspection  method  in  reducing  the  POF  or  reducing  the 
cost  of  inspections.  As  a  result,  it  is  useful  to  have  a  quantitative  estimate  how  the  probability-of- 
failure  varies  as  a  function  of  the  parameters  of  a  particular  POD  curve.  Therefore,  a 
methodology  was  derived  that  provides  a  convenient  low-cost  means  to  calculate  the  POF 
sensitivities,  that  is,  the  partial  derivative  dPJdO,  where  Pf  is  the  probability-of-failure  and  6  is 

any  parameter  of  a  POD  curve. 

1.2  Methodology  Development 

The  basic  methodology  is  summarized  here  for  a  single  inspection.  A  generalization  to  more 
inspections  can  be  found  in  [13]. 

The  probability-of-failure  is  determined  by  evaluating  the  integral 


(1) 


where  x  defines  a  vector  of  random  variables  such  as  initial  crack  size,  external  load,  crack 
growth  rate,  fracture  toughness,  fx  represents  the  joint  probability  density  function  of  the 
random  variables,  and  I(x,t)  denotes  an  indicator  function  such  that 
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(2) 


jo  g(x,?)>0  safe 
|l  g(x,t)  <  0  fail 


For  a  lifing  analysis,  g(x,t)  is  usually  defined  as  g(x,t)  =  tf  —t0,  where  tf  is  the  computed 
cycles  to  failure  and  t0  is  a  user-defined  limit.  Using  this  definition  of  g,  the  probability-of- 
failure  is  the  probability  of  the  structural  component  failing  before  the  user-defined  time  t0. 

The  probability-of-failure  with  an  inspection  can  be  determined  by  evaluating  the  integral 

pt  (0  =  I  J(x’  00  -  POD(Q,  a( y,  tx  )))/x  (x)dx  t  >  tx  (3) 

where  tx  denotes  the  time  of  inspection  1,  a  represents  the  crack  size,  and  y  represents  the 
random  variables  that  affect  the  crack  size;  y  is  a  subset  of  x .  The  sensitivity  of  the  POF  with 
respect  to  parameters  of  the  POD  curve  ( 0 )  can  be  determined  by  taking  the  derivative  of 
Equation  3  and  rearranging  to  yield 

cPf(t)  j  0  t<tx 

cQx  j  \XJ{x,t)D^{%a{ydi))CPOD1(Qva{y,ti))fK00dx  t  >  tx  <4> 

dPOD  (0  a(y,tx))  1 

where  Q  (6  ,a(y,t)  = - - — - - .  The  derivative 

?  ?  CPOD9(Q9,a( y,0) 

(PODq  (0  ,  a( y,  t] ))  /VB<;  can  be  computed  analytically  given  an  analytical  definition  of  the  POD 

curve. 

The  sensitivity  cPf ./ dO  can  be  approximated  in  the  same  manner  as  estimating  the  probability- 
of-failure  using  sampling  as 


c8, 


N 


0 

Nt 

YjKXi,OQ(QiMy  idO) 


t<tx 
t  >  tx 


(5) 


where  0[  represents  the  vector  of  parameters  of  PODv  Nx  denotes  the  samples  that  reach 
inspection  1,  i.e.,  have  not  failed  before  tx  nor  detected,  a(y;,tx)  denotes  the  crack  size  at  time  tx, 
y,  represents  the  random  variables  that  affect  crack  size  for  realization  i,  x,  represents  all  the 
random  variables  for  realization  i,  and  N  denotes  the  total  number  of  samples  in  the  simulation. 
Note,  the  summation  in  Equation  5  is  divided  by  N  instead  of  Nx  in  order  to  provide  the 
sensitivity  with  respect  to  the  total  samples  N. 
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1.3  Results  and  Discussion 


A  new  sensitivity  method  has  been  developed  that  can  provide  accurate  partial  derivatives  of  the 
probability-of-failure  with  respect  to  the  parameters  of  a  POD  curve.  The  accuracy  of  the 
sensitivities  was  verified  by  comparison  with  finite  difference  estimates  on  several  numerical 
examples.  The  new  method  computes  the  sensitivity  using  the  same  samples  used  to  calculate  the 
probability-of-failure,  hence  the  cost  is  negligible.  In  contrast,  the  finite  difference  approach  is 
extremely  laborious  in  that  a  new  probabilistic  analysis  must  be  conducted  for  each  sensitivity 
and  a  large  number  of  samples  is  required  for  accuracy. 

1.3.1  Numerical  Example 

An  academic  example  is  presented  to  demonstrate  the  methodology.  The  example  presumes  a 
defect  exists,  thus,  the  POF  has  been  normalized  to  1.  The  problem  is  fictitious  but  serves  to 
demonstrate  the  elements  of  the  method  using  an  easy  to  understand  problem.  The  problem 
consists  of  an  edge  through  crack  growing  in  a  semi-infinite  rectangular  plate  under  constant 
amplitude  loading.  The  initial  crack  size  and  material  properties  are  representative  of  titanium, 
however,  the  POD  curve  is  contrived.  Four  random  variables  were  considered:  initial  crack  size 
(lognormal),  Paris  crack  growth  constants  C  and  n  (correlated  normal)  and  fracture  toughness 
(normal).  The  problem  parameters  are  given  in  Table  1.  The  random  variable  vector  X  consists 
of  X  e  at,C,n,KIC,  and  the  Y  vector  of  random  variables  consists  of  Ye  an  C,  n . 

For  this  crack  geometry  using  the  Paris  crack  growth  law,  the  crack  size  as  a  function  of  cycles 
a,  the  critical  crack  size  ac,  the  cycles-to-failure  tf,  the  crack  size  aq  at  any  time  of  inspection 

t  ,  and  the  limit  state  g(x,t)  can  be  determined  in  closed  form  as  shown  below. 


dal  dt  =  C(AK)m 


(6) 


K, 


(7) 


ac  =  (Kic  /(1 . 12crimx  v/tt))2 


(8) 


(9) 


(10) 


g(x,t)  =  tf-t0 


(ID 
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where  t0  is  a  specified  number  of  cycles. 


The  values  used  in  the  analysis  are  shown  in  Table  1.  For  this  problem,  the  initial  crack  size, 
Paris  constants  and  the  fracture  toughness  were  random. 


Table  llnput  Parameters  for  Example  Problem 


A K 

1.12Acr 

A&(=  ^,mx  ) 

61 5  MPa 

KIC 

N[55,5.5]  MPasfm 

P Cm 

-0.9751 

Logio(C) 

A[-l  1.8,0. 157] 

m 

)V[3.81,0.146] 

ai 

Z[1 5. 1,8.48]  jam 

A  single  inspection  was  simulated  at  10,000  cycles  using  a  lognormal  POD  of  the  form 

POD(a)  =  2  +  2  ErJ 

function  of  cycles  with  and  without  inspection  is  shown  in  Figure  1 .  The  sensitivity  of  the 
probability-of-failure  with  respect  to  //,  and  <jx  as  a  function  of  cycles  is  given  in  Figures  2  and 
3,  respectively,  with  the  solid  (red)  line  denoting  the  results  using  the  equations  developed  here 
and  finite  difference  results  denoted  by  the  dotted  (black)  points.  One  million  samples  were  used 
for  the  calculation.  The  results  indicate  a  very  good  agreement  with  the  finite  difference 
solutions  verifying  the  accuracy  of  the  methodology.  Finite  difference  sensitivity  estimates 
require  a  separate  reanalysis  of  the  POF  using  a  perturbed  value  of  //j  and  ax  (2  additional 
analyses),  e.g.,  cPf  Idfi [  «  APf  / A//,  =  (Pf(jux  +  Ajux)  -P,  (//,)) /A//, .  A  large  number  of  samples 

must  be  used  to  ensure  that  the  sampling  variance  is  sufficiently  small  to  accurately  estimate  the 
difference  between  two  near-equal  probabilities.  Also,  the  accuracy  is  dependent  upon  the  value 
chosen  for  A//,.  As  a  result,  the  finite  difference  operations  are  costly  since  they  require  multiple 
analyses  and  a  large  number  of  samples.  Conversely,  using  the  methodology  presented  here,  all 
sensitivities  are  obtained  from  a  single  analysis  through  post-processing. 

The  results  for  a  single  inspection  shown  here  generalize  to  any  number  of  inspections  with 
similar  accuracy.  For  multiple  inspections,  methodologies  to  address  two  cases  were  developed 
and  demonstrated:  different  POD  curves  used  during  for  the  inspections,  and  b)  the  same  POD 
curved  used  for  all  inspections. 


In  [d\-  ju 
a/2ct  j 


with  parameters  //,=-'  9,  cr,  =0.8.  The  normalized  POF  as  a 
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Probability— of— Failure 


Inspection  at  10,000  Cycles 


Figure  1  Normalized  Probability-of-Failure  with  and  without  Inspection 
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dP  /  dp 


POF  Sensitivity  wrt  POD  Parameter 


Figure  2  Single  Inspection  Case:  Probability-of-Failure  with  respect  to  //, 
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d  P,  /  do 


0.02 


POF  Sensitivity  wrt  POD  Parameter  o1 


01  23456789  10 


Cycles  xio4 

Figure  3  Single  Inspection  Case:  Probability-of-Failure  with  respect  to  o] 
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2  Complex  Variable  Methods  for  Sensitivity  Analysis 


2.1  Introduction 

Shape  sensitivity  analysis  of  finite  element  models  is  useful  for  optimization  and  design 
modifications.  The  use  of  complex  variable  methods  for  shape  sensitivity  analysis  has  some 
potential  advantages  over  other  methods.  In  particular,  implementation  for  first  order  sensitivities 
is  straightforward  as  the  only  the  additional  requirement  is  the  analysis  of  a  finite  element  model 
containing  a  perturbation  of  the  finite  element  mesh  along  the  imaginary  axis,  that  is,  the  real 
valued  coordinates  of  the  mesh  are  unaltered.  Higher  order  sensitivities  can  be  determined  using 
analyses  with  perturbations  along  the  real  and  imaginary  axes.  Although  simple  in  concept, 
standard  finite  element  codes  (including  commercial  ones)  do  not  allow  complex  variable  nodal 
coordinates,  hence,  special  purpose  codes  were  developed  under  this  program.  The  methodology 
is  demonstrated  using  a  two  dimensional  finite  element  model  on  problems  with  known 
analytical  solutions.  It  is  found  that  the  error  in  the  sensitivities  is  primarily  defined  by  the  error 
in  the  finite  element  solution,  not  the  error  in  the  sensitivity  method. 


2.2  Methodology  Development 

Sensitivities,  otherwise  known  as  partial  derivatives,  are  easily  calculated  through  finite 
differencing  methods.  Finite  differencing  requires  that  a  function  be  evaluated  at  additional 
sample  points  along  the  real  axis  and  the  derivative  of  the  function  estimated  by  calculating  the 
relative  difference  in  the  function’s  value  divided  by  the  difference  in  the  sample  points. 

Over  the  last  twenty  years,  alternative  numerical  differentiation  techniques  have  emerged  for  use 
in  sensitivity  analysis.  Two  of  these  methods  are  complex  variable  based:  the  complex  Taylor 
series  expansion  (CTSE),  also  referred  to  as  the  complex  step  derivative  method,  and  Fourier 
differentiation  (FD).  These  methods  offer  more  accurate  and  stable  derivatives  compared  to 
finite  differencing. 

CTSE  was  first  described  by  Lyness  and  Moler  in  the  late  1960’s  [14,15].  It  reemerged  as  a  tool 
for  engineering  analysis  with  a  paper  by  Squire  and  Trapp  in  1998  [16].  Since  then  it  has  been 
used  in  a  wide  variety  of  engineering  fields  including  computational  fluid  dynamics,  dynamic 
system  optimization  and  others  [17-24].  In  all  of  these  fields,  CTSE  has  offered  a  significant 
improvement  in  accuracy  over  standard  finite  differencing  methods. 

Fourier  differentiation  was  also  developed  by  Lyness  in  the  late  60’s  and  early  70’s[  14, 15,25]. 
The  method  was  further  described  by  Henrici  and  more  recently  Bagley  [26,27].  The  method 
utilizes  additional  sample  points  in  the  complex  plane  and  an  FFT  routine  to  calculate  derivatives 
including  high  order  derivatives  with  exceptional  accuracy.  To  date,  the  method  has  not  been 
widely  used  for  the  determination  of  sensitivities  for  engineering  problems. 
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2.2.1  Numerical  Differentiation 


Numerical  differentiation  is  a  process  through  which  an  estimate  of  a  function’s  derivative  can 
be  obtained.  A  derivative  is  defined  as  the  limit  of  the  change  in  a  function’s  value  across  two 
different  points  as  the  distance  between  the  two  points  goes  to  zero. 


f(xa)  =  lim 

x— 


f(x)~  fix  „) 
x-xa 


(12) 


Finite  differencing  methods  estimate  derivatives  by  approximating  the  limit  in  Equation  12  as  a 
difference  between  a  function  evaluated  a  two  distinct  points  located  a  distance  h  apart  divided 
by  h. 


f'(x  fixo  +  h)-fixo ) 
h 


(13) 


This  distance,  h,  is  often  called  the  step  size.  When  h  is  positive,  the  method  is  referred  to  as 
forward  differencing.  When  h  is  negative  it  is  called  backwards  differencing.  When  the  forward 
difference  and  the  backwards  difference  results  are  averaged,  the  method  is  called  central 
differencing  (CD).  The  equation  for  central  differencing  is  as  follows. 


/'U)»/fc  +  *)-/(*°-/l>  (14 

J  oJ  ^ 

The  approximation  of  the  derivatives  as  a  difference  between  two  nearly-equal  numbers  leads  to 
error  due  to  the  truncation  of  terms  in  the  function’s  Taylor  series.  This  error  can  be  eliminated 
by  making  the  step  size  as  small  as  possible.  However,  as  the  step  size  gets  very  small,  a  new 
source  of  error  arises.  This  new  error  is  round-off  error  and  it  is  due  to  the  fact  that  a  computer 
cannot  accurately  calculate  a  small  difference  between  two  near-equal  numbers.  This  means  that 
for  finite  differencing  there  is  a  lower  limit  on  the  step  size  and  also  a  limit  on  the  maximum 
achievable  accuracy. 

For  the  forward  differencing  method,  all  Taylor  series  terms  above  the  first  order  term  are 
ignored.  This  means  that  the  order  of  accuracy  for  a  given  step  size  is  0(h).  Using  CD,  all  the 
even  order  terms  in  the  Taylor  series  cancel  and  the  accuracy  of  the  method  becomes  0(li2). 

Higher  order  derivatives  can  also  be  calculated  through  CD  by  using  additional  sample  points. 
The  formula  for  the  second  derivative  is. 


f2\x0) 


f(xa  +  h)-2f(x0)  +  fjxa  -  h) 
h2 


(15) 
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where  the  superscript  (2)  denotes  the  second  derivative. 

One  of  the  problems  with  CD  is  that  the  calculation  of  higher  order  derivatives  requires  more 
sample  points  and  more  difference  operations.  Each  additional  difference  operation  results  in  an 
increase  in  the  round-off  error,  which  further  restricts  the  lower  limit  of  h.  This  means  that  CD  is 
not  a  good  choice  for  the  calculation  of  higher  order  derivatives. 

CTSE  is  another  numerical  differentiation  method  similar  in  concept  to  finite  differencing.  CTSE 
uses  the  orthogonality  of  the  real  and  imaginary  axes  of  the  complex  plane  to  calculate 
derivatives  with  fewer  difference  operations  and  in  turn  less  round-off  error  when  compared  to 
CD.  Similar  to  finite  differencing,  CTSE  requires  the  difference  of  two  analyses  but  with  a  small 
perturbation  along  the  imaginary  axis.  That  is,  variable  X  =  x0  is  perturbed  to  X  =  x0  +ih ,  where 
i  denotes  an  imaginary  number  and  h  denotes  the  step  size.  The  formulae  for  the  derivatives  can 
be  derived  from  the  Taylor  series  representation  of  the  function  evaluated  at  the  complex  sample 
point. 

/Oo  + ih)  =  /(*o)  +  fa\xo)l-^+f(2\xo)^-  +  /(3,(*0)^  +  L  <16) 

where  /(1)  denotes  the  first  derivative,  /(2)  the  second,  etc.  Taking  the  imaginary  part  of  both 
sides  of  Equation  16  and  solving  for  the  first  derivative  will  result  in  an  approximation  with 
accuracy  O(lr). 


wiw  .  _  /(x0+z7z)-/(x0)  ^  I  m(/(x„ +//?)) 
J  ^  °’~  h  ~  h 


(17) 


It  is  noted  that  no  difference  operation  is  needed  for  the  first  derivative.  This  means  that  the  step 
size  can  be  made  arbitrarily  small  with  no  concern  about  increasing  round-off  error.  Taking  the 
real  part  of  Equation  16,  the  formula  for  the  second  derivative  with  error  0(h)  can  be  derived. 

f»(x  )aZ(/(v)-R</(v+ft)))  (18) 

h 

It  is  noted  that  the  second  derivative  contains  a  difference  operation  meaning  that  round-off  error 
will  be  a  problem  if  h  is  set  too  small.  By  using  more  sample  points  along  the  imaginary  axis  it  is 
possible  to  solve  Equation  16  to  obtain  the  higher  order  derivatives  [28]. 

Higher  order  derivatives  can  be  computed  using  complex  variable  sensitivity  method  such  as 
Fourier  Differentiation  (FD)  [27].  The  heart  of  Fourier  differentiation  is  making  a  real  valued 
function  become  a  periodic  complex  function.  This  is  accomplished  by  adding  a  periodic, 
oscillatory,  complex  component  to  each  of  the  function's  real  independent  variables.  The 
resulting  periodic  function  now  has  a  Taylor  series  representation  that  takes  on  the  properties  of 
a  periodic  Fourier  series.  Fast  Fourier  transform  techniques  are  used  to  determine  the  coefficients 
of  this  series.  The  resulting  coefficients  contain  the  function's  derivatives. 
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A  central  feature  of  Fourier  differentiation  is  to  make  the  individual  terms  in  the  Taylor  series 
oscillate  at  different  frequencies.  To  achieve  these  oscillations,  the  perturbation  to  the  variable  x 
is  taken  to  be  a  complex  number  as  Ax  =  ce  '  ,  where  c  is  a  sampling  radius.  This  expression 
extracts  from  the  series  the  nth  coefficient  that  contains  the  nth derivative  of  the  function, 

1  f  ft  ,  ^  w  fn\xo)cn  ri«, 

—  J  f{xa  +ce  )e  d0  = - - -  (19) 

2/zr  n\ 


If  the  function  of  interest  is  evaluated  at  N  sample  points  along  a  circular  contour  c  in  the 
complex  plane  centered  on  the  initial  point,  a  vector  of  the  sampled  data  can  be  run  through  an 
FFT  routine  and  the  output  will  be  the  first  N  terms  in  the  function’s  Taylor  series.  The  nth  order 
derivative  of  the  function  can  then  be  calculated  from  the  Taylor  series  coefficients  by  using  the 
following  relationship. 


a  in 


(20) 


where  a„  is  the  nth  Taylor  series  coefficient.  For  more  information  on  Fourier  differentiation  see 
ref.  [27], 


2.3  Results  and  Discussion 

2.3.1  Numerical  Example:  Thick  Walled  Cylinder 

Two  problems  with  analytical  solutions  were  examined  in  order  to  test  the  accuracy  of  the  three 
numerical  differentiation  techniques;  a  thick  walled  cylinder  under  uniform  boundary  pressure 
and  a  disk  under  diametrical  compression.  The  equations  that  govern  the  stress  through  the 
thickness  of  the  cylinder  are  given  in  Equation  21  [29], 


2  2,  2 

_  n  r2  pi  r2p 

^  r  2  2  2  2  2 

ri  ~n  r  ri  ~n 

2  2,  2 

n  r2  p  1  r2  p 


<Jn=- 


2  2  2 

r2  -rx  r 


r2  -rx 


(21) 


where,  r/  is  the  inner  radius,  r2  is  the  outer  radius,  and  p  is  the  boundary  pressure.  For  this 
example,  the  inner  radius  of  the  cylinder  is  0.5  m,  the  outer  radius  is  1  m  and  the  boundary 
pressure  is  10  kPa.  The  stress  equations  given  in  Equation  21  can  be  differentiated  with  respect 
to  the  inner  radius  to  generate  the  sensitivities  of  the  stresses.  The  first  and  second  order 
sensitivities  of  the  stresses  with  respect  to  the  inner  radius  appear  in  Equations  22  and  23. 
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for 

drx 


=  2 


n  riP 

/  2  2x2 

(h  -n ) 


\ 

l 

J 


fog 

r\r2  P 

(  2  3 

^  +1 

drx 

(  2  2x2 

Lte  ) 

^  )\ 

*°'=2 

(3?i  +  r22)r22p 

(  2  \ 

\-l 

3\ 

/  2  2\  3 

(r2  -n ) 

l r  )_ 

dt\~ 


(22) 


(23) 


This  problem  was  solved  using  four  different  meshes  in  order  to  examine  the  convergence  of  the 
error  in  the  solution.  The  coarsest  mesh  contained  888  elements  and  1868  nodes;  the  next  mesh 
contained  1792  elements  and  3716  nodes;  the  third  mesh  contained  6184  elements  and  12,608 
nodes;  and  the  finest  mesh  consisted  of  25,124  elements  and  50,724  nodes.  The  solutions  for  the 
radial  and  tangential  stresses  as  calculated  using  the  6184  element  mesh  are  shown  in  Figure  4. 


Figure  4.  The  Numerical  Solution  of  the  Radial  and  Tangential  Stresses  for  Example  1.  A) 
The  FEM  Solution  for  the  Radial  Stresses  B)  The  FEM  Solution  for  the  Tangential  Stresses 

The  following  norm  was  selected  in  order  to  compare  the  error  in  the  four  different  mesh  cases. 


I  .I  \<Janaiyticai  O 'numerical ) 

\error\\  = - -r - j - 

mean(\(7analytica\) 


(24) 
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Table  2  shows  the  norm  of  the  error  in  both  the  radial  and  tangential  stress  solutions  for  each 
mesh  case.  This  data  is  shown  graphically  in  Figure  5.  It  is  seen  that  each  successive  mesh 
iteration  reduces  the  error  by  approximately  half  an  order  of  magnitude. 

Table  2.  The  Norm  of  the  Error  in  the  Stress  Solutions  for  a  Thick  Walled  Cylinder 


Number  of  Elements 

Radial  Stress 

Tangential  Stress 

888 

5.251E-3 

3.1852E-3 

1792 

2.626E-3 

1.7412E-3 

6184 

7.574E-4 

5.3097E-4 

25124 

1.860E-4 

1.3248E-4 

Figure  5.  Convergence  of  the  Error  in  the  Radial  and  Tangential  Stress  Models  for  Thick 
Walled  Cylinder.  A.  The  Norm  of  the  Error  in  the  Radial  Stress,  B.  The  Norm  of  the  Error 

in  the  Tangential  Stress 

The  amount  of  computational  time  (wall  time)  needed  to  solve  each  finite  element  model  is 
shown  in  Table  3.  Given  the  amount  of  computational  time  required  for  the  25,124  element  case 
and  the  fact  that  the  complex  sensitivity  solutions  will  require  three  times  more  computation  than 
the  real  valued  case,  the  6184  element  case  was  used  to  generate  the  first,  second  and  third  order 
sensitivities.  For  each  sensitivity  method  the  step  size  or  sampling  radius  was  0.001,  which  is 
approximately  l/30th  of  the  average  element  edge  length.  CTSE  and  CD  were  both  performed 
using  as  few  sample  points  as  possible,  and  FD  was  performed  using  6  sample  points. 
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Table  3.  The  Computational  Time  Required  to  Solve  each  Model  for  a  Thick  Walled 
_ Cylinder _ 


Number  of  Elements 

Time  For  Solution  (s) 

888 

8.49 

1792 

19.48 

6184 

130.76 

25124 

2799.77 

The  pointwise  error  in  the  first  and  second  order  sensitivities  is  shown  in  Figure  6  and  7, 
respectively  using  the  following  formula, 


analytical  numerical  , 

error  = - : — j - i —  (25) 

maX(|  °analy,ica\) 

This  formula  was  chosen  so  as  to  minimize  the  influence  of  large  errors  at  locations  where  the 
solution  is  near  zero  such  as  at  the  inner  radius  for  the  radial  stress. 

The  norm  of  the  error  in  the  first,  second,  and  third  order  sensitivities  appears  in  Table  4.  The 
error  in  the  first  order  sensitivities  of  the  radial  stress  over  the  entire  domain  appear  in  Figure  6, 
and  the  error  in  the  second  order  sensitivities  of  the  radial  stress  appear  in  Figure  7.  These  figures 
show  only  very  slight  differences  between  the  three  methods.  It  is  also  seen  that  along  the  inner 
circumference  of  the  cylinder  the  error  is  large.  This  is  due  to  the  fact  that  the  sensitivity  cannot 
be  accurately  calculated  on  the  inner  radius  because  the  boundary  conditions  require  the  solution 
to  be  fixed  at  the  inner  surface. 


Table  4.  The  Norm  of  the  Error  in  the  Sensitivity  of  the  Stress  to  the  Inner  Radius  for  a 
_ _ Thick  Walled  Cylinder _ 


Method 

Radial  Stress 

Tangential  Stress 

Order 

1st 

2nd 

3rd 

1st 

2nd 

3rd 

CD 

4.7265E-2 

8.1689E-2 

2.4292E-2 

8.4286E-2 

3.2213E-1 

1.4765E-2 

CTSE 

4.7268E-2 

8.1766E-2 

2.4298E-2 

8.2488E-2 

3.1587E-1 

1.4769E-2 

FD 

4.7268E-2 

8.1671E-2 

2.4295E-2 

8.7353E-2 

3.3258E-1 

1.4769E-2 
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Figure  6.  The  Error  in  the  First  Order  Sensitivity  of  the  Radial  Stress  for  Example  1.  A) 
Error  in  CD,  B)  error  in  CTSE,  C)  error  in  FD. 
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Figure  7.  The  Error  in  the  Second  Order  Sensitivity  of  the  Radial  Stress  to  the  Inner 
Radius  for  Example  1.  A)  Error  in  CD,  B)  error  in  CTSE,  C)  error  in  FD. 

The  norm  of  the  error  in  each  case  is  mostly  independent  of  the  method  selected.  This  is 
especially  true  for  the  first  order  sensitivities.  The  fact  that  the  error  is  similar  between  all  three 
methods  points  to  the  fact  that  the  error  in  the  solution  is  dominating  errors  arising  from  the 
differentiation  methods  themselves.  This  is  seen  by  looking  at  the  first  and  second  order 
sensitivities  of  the  radial  stress  as  a  function  of  the  number  of  elements  shown  in  Table  5.  It  is 
quickly  seen  that  each  additional  mesh  refinement  increases  the  accuracy  of  the  method.  This 
reduction  in  the  errors  of  the  sensitivities  is  similar  to  the  reduction  in  the  error  of  the  solution 
due  to  further  mesh  refinement  seen  in  Table  2. 
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Table  5.  The  Norm  of  the  Error  in  the  First  Order  Sensitivity  of  the  Radial  Stress  as  a 
Function  of  the  Number  of  Elements  for  a  Thick  Walled  Cylinder 


Number  of 
Elements 

CD 

CTSE 

FD 

Order 

1st 

2nd 

1st 

2nd 

1st 

2nd 

888 

1.27E-1 

2.08E-1 

1.27E-1 

2.08E-1 

1.27E-1 

2.08E-1 

1792 

8.97E-2 

1.53E-1 

8.99E-2 

1.52E-1 

8.99E-2 

1.52E-1 

6184 

4.72E-2 

8.16E-2 

4.72E-2 

8.17E-2 

4.72E-2 

8.16E-2 

25124 

2.29E-2 

4.00E-2 

2.29E-2 

4.00E-2 

2.29E-2 

3.98E-2 

The  error  in  the  first  and  second  order  sensitivities  of  the  radial  stress  calculated  using  three 
different  step  sizes  is  shown  in  Table  6.  It  is  seen  that  changing  the  step  size  does  not  have  much 
effect  on  the  accuracy  of  the  sensitivity.  This  is  a  further  indicator  that  the  accuracy  of  the 
solution  is  limiting  the  accuracy  of  the  sensitivities,  not  the  accuracy  of  the  numerical 
differentiation  methods.  One  exception  is  the  second  order  sensitivity  at  the  smallest  step  size, 
0.0001  or  l/300th  of  the  average  element  edge  length.  At  this  step  size  each  method  produces 
sensitivities  that  are  less  accurate  than  those  calculated  with  a  larger  step  size.  This  indicates  that 
the  machine  round-off  error  associated  with  this  step  size  may  be  similar  in  magnitude  to  the 
error  due  to  the  solution. 

Table  6.  The  Norm  of  the  Error  in  the  First  and  Second  Order  Sensitivities  of  the  Radial 
Stress  as  a  Function  of  Step  Size  for  a  Thick  Walled  Cylinder 


Step  Size 

CD 

CTSE 

FD 

Order 

1st 

2nd 

1st 

2nd 

1st 

2nd 

.0001 

4.7268E-2 

8.2959E-2 

4.7268E-2 

1.0500E-1 

4.7268E-2 

9.2942E-2 

.001 

4.7265E-2 

8.1689E-2 

4.7268E-2 

8.1766E-2 

4.7268E-2 

8.1671E-2 

.01 

4.7091E-2 

8.2920E-2 

4.7249E-2 

8.1197E-2 

4.7271E-2 

8.1748E-2 

2.3.2  Numerical  Example:  Disc  in  Diametrical  Compression 

One  of  the  classical  tests  in  material  analysis  is  the  disc  in  diametrical  compression  [20].  In  this 
test  a  circular  disc  is  loaded  in  compression  along  its  y-axis.  The  load  is  modeled  as  a  point  load. 
This  loading  and  geometry  generates  a  very  nice  uniform  tensile  stress  along  the  x-axis  of  the 
specimen.  It  is  thus  useful  in  examining  the  tensile  properties  of  a  material  without  actually 
loading  the  specimen  in  tension.  This  test  is  also  known  as  the  indirect  tensile  test  or  the  Brazil 
nut  test. 

The  diametrical  compression  test  has  an  analytical  solution  that  can  be  derived  through  simple 
superposition.  The  solution  of  the  stresses  is  given  in  Equation  26  [29]. 
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(26) 


In  these  equations,  P  is  the  magnitude  of  the  point  load,  R  is  the  radius  of  the  disc,  and  x  and  y 
specify  the  location  at  which  the  stress  is  calculated,  with  the  point  (0,0)  located  at  the  center  of 
the  disc.  The  analytic  solutions  of  Equation  26  can  be  differentiated  to  yield  the  sensitivities  with 
respect  to  the  radius  of  the  disc.  The  equations  for  the  first  two  sensitivities  of  the  normal  stress 
in  the  x-direction  with  respect  to  the  radius  are 


d°x 

dR 


2 P 
n 


(3 R1  -  6Ry  -x2  +  3 y2)x2  (3R2 +  6Ry  -  x2 +  3y~)x2 
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12  (R -y)(R2 -2  Ry -x2  +y2)x2 
-2 P  (R2  -2Ry  +  x2  +y2)4 

n  12(7?  +  y)(R2  +2Ry -x2  +y2)x2  1 

_  (R2  +2Rv  +  x2  +y2)4  R2 


(27) 


The  closed-form  solutions  for  the  sensitivities  make  this  problem  another  excellent  choice  for 
exploring  the  use  of  the  complex  variable  sensitivity  methods. 

The  diametrical  compression  test  model  was  solved  using  three  different  meshes,  with  a  coarse 
mesh  consisting  of  1148  elements  and  2357  nodes;  a  moderately  refined  mesh  of  2502  elements 
and  5093  nodes;  and  a  fine  mesh  with  8,374  elements  and  16,909  nodes.  The  solution  for  the 
stresses  as  calculated  using  the  fine  mesh  appears  in  Figure  8.  It  should  be  noted  that  for  each 
figure  in  this  example,  no  solution  is  plotted  for  the  elements  that  share  the  node  where  the  load 
is  applied.  This  is  due  to  the  fact  that  the  stress  on  this  node  would  be  infinite. 
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Figure  8.  The  Finite  Element  Solution  for  the  Stresses  in  a  Disc  in  Diametrical 
Compression.  A)  The  Numerical  Solution  for  the  Stresses  in  the  x-direction  B)  The 
Numerical  Solution  for  the  Stresses  in  the  y-direction  C)  The  Numerical  Solution  for  the 

Shear  Stress 


The  error  norm  used  for  this  example  is 


error  = 


® analytical  ^ numerical 


analytical 


(28) 


Table  7  shows  the  error  for  the  three  stresses  for  the  three  different  mesh  sizes  using  the  error 
norm  give  in  Equation  18.  As  before,  it  is  seen  that  each  successive  mesh  refinement  results  in  a 
significant  reduction  in  the  error  norm.  This  is  shown  visually  in  Figure  9.  The  red  areas  are 
regions  of  relatively  higher  error.  As  the  number  of  elements  increases,  the  total  size  of  the  red 
regions  decreases  significantly  as  the  mesh  is  further  refined.  The  computational  time  required  to 
generate  one  solution  for  the  three  different  mesh  sizes  is  shown  in  Table  8. 
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Table  7.  The  Norm  of  the  Error  in  the  Stress  Solutions  for  a  Disk  under  Diametrical 
_ _ Compression _ _ 


Number  of 
Elements 

Norm  of  Error  in  Stress 
in  X 

Norm  of  Error  in 
Stress  in  Y 

Norm  of  Error  in 
Shear  Stress 

1148 

1.261E-1 

4.379E-2 

9.717E-2 

2502 

8.633E-2 

3.234E-2 

6.122E-2 

8374 

4.578E-2 

2.085E-2 

3.953E-2 
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Figure  9.  The  Error  in  ax  for  Three  Different  Meshes  for  a  Disc  in  Diametrical 
Compression.  A)  The  error  for  the  mesh  with  1,148  elements  and  2,357  nodes,  B)  The  error 
for  the  mesh  with  2,502  elements  and  5,093  nodes,  C)  The  error  for  the  mesh  with  8,374 

elements  and  16,909  nodes. 
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Table  8.  The  Computational  Time  Required  to  Solve  each  Model  for  a  Disk  under 
_ Diametrical  Compression _ 


Number  of  Elements 

Solution  Time  (s) 

1148 

10.41 

2502 

27.82 

8374 

200.00 

The  errors  in  the  sensitivities  of  the  normal  stress  in  the  x-direction,  calculated  by  each  of  the 
three  methods  are  plotted  in  Figure  10  (first  order)  and  Figure  1 1  (second  order)  for  the  fine 
mesh  case.  It  is  seen  that  there  is  again  relatively  higher  error  along  the  circumference  of  the  disc 
due  to  boundary  conditions.  It  is  also  noted  that  a  few  lines  of  relatively  high  error  form  inside 
the  discs.  These  lines  represent  regions  where  the  analytical  sensitivity  is  zero  or  near  zero.  Since 
the  analytical  sensitivity  appears  in  the  denominator  of  the  error  formula  given  in  Equation  28, 
the  error  becomes  very  large  when  the  analytical  sensitivity  tends  towards  zero. 


1st  Derivative  of  Sigma  x  CD 


1st  Derivative  of  Sigma  CTSE 


1st  Derivative  of  Sigma  :  FD 
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Figure  10.  The  Error  in  the  First  Order  Sensitivity  for  a  Disc  in  Diametrical  Compression. 
A)  Solution  calculated  by  CD,  B)  Solution  calculated  by  CTSE,  C)  Solution  calculated  by 

FD 
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2nd  Derivative  of  Sigma  :  CO  2nd  Derivative  of  Sigma  x  CTSE 


2nd  Derivative  of  Sigma  x  FD 


Figure  11.  The  Error  in  the  Second  Order  Sensitivity  for  a  Disc  in  Diametrical 
Compression.  A)  Solution  calculated  by  CD,  B)  Solution  calculated  by  CTSE,  C)  Solution 

calculated  by  FD 

The  norm  of  the  error  in  the  sensitivities  is  shown  in  Table  9.  For  this  example,  as  compared  to 
the  thick  walled  cylinder,  there  is  not  a  large  dependence  of  the  error  norm  on  the  number  of 
elements.  Also,  there  is  not  much  difference  between  the  methods  themselves.  This  points  to  the 
fact  that  the  error  in  the  derivatives  is  not  due  to  the  truncation  error  of  the  derivative  methods, 
since  the  truncation  error  of  CTSE  and  CD  should  be  of  the  order  0(h  )  while  the  truncation 
error  of  FD  should  be  0(h6).  The  lack  of  dependence  on  the  choice  of  method  is  further  shown  in 
Table  10  where  the  norm  of  the  errors  is  shown  for  the  2502  element  mesh  for  three  different 
step  sizes,  or  sampling  radii.  Very  little  variation  in  the  error  is  seen  as  a  function  of  the  step  size, 
which  is  not  the  behavior  of  truncation  error. 
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Table  9.  The  Norm  of  the  Error  in  the  First  Order  Sensitivity  of  ax  as  a  Function  of  the 
Number  of  Elements  for  a  Disk  under  Diametrical  Compression 


Number  of  Elements 

CD 

CTSE 

FD 

Order 

1st 

2nd 

1st 

2nd 

1st 

2nd 

1148 

0.4681 

0.6378 

0.4681 

0.6379 

0.4681 

0.6378 

2502 

0.4164 

0.5673 

0.4164 

0.5675 

0.4164 

0.5674 

8374 

0.4267 

0.6572 

0.4270 

0.6581 

0.4270 

0.6577 

Table  10.  The  Norm  of  the  Error  in  the  First  and  Second  Order  Sensitivities  of  the  Radial 
Stress  as  a  Function  of  Step  Size  for  a  Disk  under  Diametrical  Compression 


Step  Size 

CD 

CTSE 

FD 

Order 

1st 

2nd 

1st 

2nd 

1st 

2nd 

.0001 

0.4164 

0.5679 

0.4164 

0.5726 

0.4164 

0.5688 

.001 

0.4164 

0.5673 

0.4164 

0.5675 

0.4164 

0.5674 

.01 

0.4143 

0.5609 

0.4166 

0.5707 

0.4164 

0.5674 

2.3.3  Numerical  Example:  Sensitivity  of  stress  intensity  factor  with  respect  to  crack 
length 

In  this  example,  CTSE  was  applied  to  several  fracture  mechanics  problems  to  determine  the 
sensitivity  of  the  stress  intensity  factor  with  respect  to  crack  length.  Determining  the  stress 
intensity  factor  ( K)  at  a  crack  tip  and  its  sensitivity  to  changes  in  crack  length  (a)  are  important 
problems  in  fracture  mechanics.  A  common  method  to  calculate  the  stress  intensity  factor  from 
finite  element  method  results  is  to  numerically  compute  a  /-integral  using  the  domain  integral 
method  and  a  truncated-pyramid  virtual  displacement  (Q)  function.  As  shown  below,  the  first 
order  derivative  of  the  J  integral  (or  K)  with  respect  to  crack  length  can  be  determined  using 
CTSE  by  adding  an  imaginary  step  ( h )  to  the  elements’  nodal  coordinates  and  applying  Equation 
17.  This  methodology  achieves  high  accuracy  over  a  wide  range  of  imaginary  step  sizes  for  the 
two-dimensional  (2D)  problems  tested  to  date. 

Four  well-known  two-dimensional  fracture  mechanics  problems,  the  single-  and  double-edge 
crack,  central  crack,  and  an  infinite  array  of  collinear  cracks,  can  all  be  modeled  using  the  same 
geometry,  material  properties,  and  remote  loading,  by  varying  the  boundary  conditions  along 
symmetry  planes  as  shown  in  Figure  12. 
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a 


- Free  Surface  ,  , 

c  t  ni  Applied  Pressure 

-  Symmetry  Plane _  '  ' _ 

Figure  12.  The  Four  Problems  Modeled  for  Example  3,  A)  Single  Edge  Crack  (SENT),  B) 
Double  Edge  Crack  (DENT),  C)  Center  Crack  (CCT),  D)  Infinite  Array  of  Collinear 
Cracks.  The  darker  shaded  area  indicates  the  common  geometry. 
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The  mode  I  reference  solutions  for  these  problems  are  shown  in  Table  11.  The  infinite  array  of 
collinear  cracks  is  particularly  attractive  as  it  has  closed-form  solution  for  K.  The  other  K 
solutions  are  approximate  fits  to  the  numerical  solutions. 


Table  11  Reference  Formulas  for  the  Stress  Intensity  Factor  (K) 


Geometry 

^I,ref  ~  O'remote'JaY 

Single  Edge 
Crack[30] 
(SENT) 

„  7WL  , 

2tan - (  r  y3 

<7  V  2W  0.752+2.02  a  +0.37  1  sin  m 

remote  ml  w  \  2W  J 

cos  v  y 

2W 

Double  Edge 
Crack[31] 
(DENT) 

i —  )W {  m  ^ ,  .  2m 

Vremote-V  ™  J  tan  +  0.1  Sin 

\[  m\  W  W  y 

Central  Crack[32] 
(CCT) 

) - '  (  V 

, —  m  a 

-y\ -rn  roe  1  H  1  H 

f  V) 
a 

^remo,e^7ra']jSeC2w\^  0‘02\^J  +  0’06 

kWJ  J 

Infinite  Array[33] 

I —  2  W  m 

(7remote^mi\  tm 

V  m  2  W 

All  of  these  problems  assume  an  infinite  longitudinal  length  (L),  i.e.,  the  height  in  Figure  12,  but 
a  model  with  a  finite  geometry  produces  a  good  approximation  of  the  stress  fields  when  the  crack 
length  is  much  less  than  height  of  the  model. 

The  sensitivities  of  K  to  crack  length  can  be  found  by  differentiating  these  formulas  with  respect 
to  a  as  shown  in  Table  12. 
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Table  12  The  First  Partial  Derivative  of  Each  Reference  Formula  with  Respect  to  Crack 

Length 


K'  =  EK  lrcf  /  di 

Single  Edge  Crack 
(SENT) 

Too  lengthy  to  be  included  here 

Double  Edge  Crack 
(DENT) 

2  m  _  2  tta 

rrz  COS - b5sec~ - 

vlO  w  W 

-  //  f  t  - 

-t  remote  1  /  x 

10  urf  ■  2MI  Ml) 

W  sin - h  10  tan  — 

V  l  w  wj 

Central  Crack 
(CCT) 

Cr remote  1  ^  "qq  ^ 

800 W5  V  a  “  2W 

^4001T5  +200^W4<atan  m  50Vk3a2 

2W 

-57rW2a3  tanJ^  +  2161Tfl4  +12^<25  tan^^ 
l  2W  2  W. 

Infinite  Array 

2  m 

rg  sec 

2W 

-  TTC7  - 

A  remote  1 

4  TT  T  71(1 

.  W  tan - 

V  2  IT 

In  these  problems,  the  stresses  around  the  crack  tip  are  pure  mode  I  so  the  stress  intensity  factor 
may  be  determined  from  the  /-integral  values  calculated  using  the  finite  element  results  as 
K,  =  a [JE'  where  E' =  E  for  plane  stress  and  E'  =  E  /(I  -  v2)  for  plane  strain. 

These  four  related  problems  use  the  same  finite  element  mesh  of  eight-noded  bi-quadratic 
quadrilaterals  with  degenerate  quarter-point  nodes  around  the  crack  tip.  The  material  properties 
and  loading  conditions  are  also  the  same  and  the  problems  differ  only  in  the  boundary  conditions 
applied.  The  numerical  values  used  for  the  problems  are  shown  in  Table  13. 

Table  13  Parameters  of  the  Plane  Strain  Problems  Investigated 


Remote  Load 

^remote 

100xl0J 

Elasticity 

E 

10xl0y 

Poisson’s  Ratio 

Y 

0.3 

Height 

L 

2 

Thickness 

B 

1 

Width 

W 

1 

Crack  Length 

A 

0.100 

In  order  to  determine  the  derivative  of  K  with  respect  to  crack  size,  a  complex  variable  finite 
code  (CFEM)  was  written  for  2D  linear  elasticity  analysis  with  a  J  integral  domain  integration 
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calculation.  The  J  integral  results  from  CFEM  were  verified  against  results  from  Abaqus  using 
identical  finite  element  meshes. 

The  recommended  mesh  to  use  around  a  crack  tip  in  order  to  obtain  accurate  J  integral 
calculations  is  called  a  “spider  web  mesh”  of  which  a  schematic  is  shown  below.  The  mesh 
consists  of  a  number  of  concentric  rings  of  elements  around  the  crack  tip.  The  mesh  in  Figure  13 
consists  of  6  concentric  rings  of  elements.  The  heavy  black  line  with  gray  interior  indicates  that 
the  J  integral  was  computed  using  this  integration  domain,  i.e.,  5  rings  of  elements.  The  elements 
in  ring  1  are  collapsed  quadrilaterals  elements  with  quarter  point  locations  for  mid-side  nodes. 
The  parameters  used  in  the  following  numerical  examples  are  given  in  Table  14. 


Figure  13  Schematic  of  "Spider  Web"  Finite  Element  Mesh  Used  for  J  Integral 

Calculations 


Table  14  Selected  Finite  Element  Model  Characteristics 


Approximate  Domain  Qi  Radius 

1.5xl0'6 

Largest  Domain  (Q36)  Radius 

0.15a 

Radial  Elements 

36 

Circumferential  Elements 

36 

Number  of  Elements 

12,482 

Number  of  Nodes 

37,957 

Degrees  of  Freedom 

75,914 

Although  complex  number  operations  typically  require  three  times  the  computational  effort  of 
their  real-valued  equivalent,  for  many  problems  only  a  limited  number  of  scalar  values  have  an 
imaginary  step  applied.  When  estimating  shape  sensitivities,  the  number  of  complex  operations 
needed  is  approximately  proportional  to  the  fraction  of  the  node  ordinates  affected  by  the 
imaginary  displacement,  so  that  the  total  computational  effort  is  less  than  it  would  be  if  all 
operations  needed  to  be  complex.  Typical  timing  estimates  are  shown  in  Table  15. 
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Table  15  Typical  Computation  Times.  Note  These  Values  Should  Not  Be  Compared  To 
Those  In  Table  4  Since  This  Example  Uses  A  Different  CFEM  Implementation 


Imaginary 

Displacement 

Meethod 

Time  to  Build 
Stiffness 
Matrix  (s) 

Time  to  Solve  for 
System 

Displacement  (s) 

Time  to  Calculate  J 
for  all  Integration 
Domains  (s) 

Total 
Solution 
Time  (s) 

None  (Pure  Real) 

27.9 

2.4 

1.6 

31.9 

Tip  Node 

28.0 

4.9 

1.7 

34.6 

Domain  1 

28.0 

4.9 

1.6 

34.5 

Crack  Face 

28.0 

4.9 

1.6 

34.5 

Domain  4 

28.0 

4.9 

1.8 

34.7 

When  a  2D  finite  element  mesh  represents  a  solid  body  with  a  straight  crack,  there  are  many 
ways  to  perturb  the  nodes  to  increase  the  crack  length  by  a  small  imaginary  amount  ih.  One 
simple  method  is  to  move  only  the  node  at  the  crack  tip  a  small  amount  in  the  direction  of  crack 
extension.  Another  method  is  to  move  all  the  nodes  along  the  crack  face  in  the  direction  of  crack 
extension  proportional  to  the  distance  from  the  crack  mouth  so  that  the  tip  node  moves  ih  and  the 
nodes  at  the  mouth  remain  unperturbed.  These  methods  have  the  disadvantage  of  changing  the 
relative  locations  of  some  of  the  effected  elements’  mid-side  nodes,  i.e.,  the  mid-side  nodes  on 
the  crack  face,  so  they  are  no  longer  precisely  at  mid-  or  quarter-points,  and  therefore  represent  a 
change  in  both  crack  length  and  a  change  to  the  elements’  shape  functions,  resulting  in  less 
accurate  results.  Although  the  effect  is  not  large,  it  does  affect  the  accuracy  of  the  computed 
derivative  of  the  J  integral.  A  schematic  of  this  effect  is  shown  below. 


Figure  14  Schematic  of  the  Effect  of  Perturbing  only  the  Nodes  Along  the  Crack  Face 


A  solution  to  this  problem  is  to  perturb  all  the  elements  within  a  user-defined  group  around  the 
crack  tip  as  a  rigid  body  in  the  imaginary  domain.  For  example,  in  Figure  15,  all  the  nodes  along 
the  black  line  (ring  5)  remain  unperturbed,  whereas,  the  nodes  within  ring  4  are  all  perturbed  as  a 
rigid  body.  In  this  case,  we  call  the  perturbation  method  Q4  to  indicate  the  all  the  nodes  within 
ring  4  are  perturbed  as  a  rigid  body.  Obviously,  any  ring  can  be  chosen  for  the  perturbation.  As 
shown  in  figures  below,  the  highest  accuracy  appears  to  be  obtained  when  using  Q4  for  the 
imaginary  perturbation,  then  using  the  J  integral  results  obtained  from  a  ring  of  elements  larger 
than  Q4. 
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Figure  15  Schematic  of  Rigid  Body  Perturbation  of  Nodes  Around  the  Crack  Tip 

Since  the  effect  of  the  complex  step  h  is  relative  to  the  size  of  the  elements  having  their  nodes 
perturbed,  the  complex  step  h  is  normalized  by  h\,  with  hi  defined  as  the  smallest  real-valued 
distance  by  which  any  comer  node  can  be  moved  in  the  displacement  direction  without 
overtaking  another  node. 

Table  16  compares  the  mode  I  stress  intensity  factor  and  the  stress  intensity  factor  derivative 
results  for  the  infinite  array  problem  from  CFEM  using  the  J  integral  ( KCFEM ),  computed  using 
the  J  integral  with  Abaqus  (KAb  ),  from  the  reference  solution  (Kref),  and  the  first  order 
derivatives  obtained  using  CTSE  within  CFEM  (K'CFEM)  and  the  solution  in  Table  12  (K’ref).  Four 

different  methods  were  used  to  implement  the  perturbation  of  the  crack  length  in  the  imaginary 
domain.  The  results  reported  are  obtained  from  the  36th  ring  of  circumferential  elements. 
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Table  16  Accuracy  of  K  and  K'  for  Four  Methods  of  Imaginary  Displacement  for  the 
Infinite  Array  Problem,  J-Integration  Domain  36,  h/hi=  1CT6 


^CFEM^Abq 

^CFEM^ref 

K'cfem  !  K'ref 

Single  Edge  Crack 

Tip  Node 

0.983060 

0.977694 

0.990846 

Domain  1 

0.983060 

0.977694 

0.990833 

Crack  Face 

0.983060 

0.977694 

0.991554 

Domain  4 

0.983060 

0.977694 

0.989891 

Double  Ed 

ge  Crack 

Tip  Node 

0.999718 

1.010374 

0.962984 

Domain  1 

0.999718 

1.010374 

0.962971 

Crack  Face 

0.999718 

1.010374 

0.963675 

Domain  4 

0.999718 

1.010374 

0.962056 

Central  Crack 

Tip  Node 

0.999510 

0.999613 

1.001053 

Domain  1 

0.999510 

0.999613 

1.001039 

Crack  Face 

0.999510 

0.999613 

1.001776 

Domain  4 

0.999510 

0.999613 

1.000088 

Infinite  Array 

Tip  Node 

0.999591 

0.999623 

1.000761 

Domain  1 

0.999591 

0.999623 

1.000748 

Crack  Face 

0.999591 

0.999623 

1.001485 

Domain  4 

0.999591 

0.999623 

0.999797 

Figure  16  shows  the  results  for  K  as  a  function  of  the  J  integration  domain.  The  results  indicate 
that  the  most  accurate  solution  is  obtained  from  the  J  integral  domain  4,  e.g.,  4th  ring  of  elements. 
This  result  appears  to  be  independent  of  the  complex  perturbation  method. 
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J  Integration  Domain 

Figure  16.  CFEM  and  ABAQUS  Infinite  Array  Stress  Intensity  Factor  (K)  Accuracy  vs.  J 

integral  domain 

The  results  in  Table  16  and  Figure  16  indicate  close  agreement  in  K  between  CFEM  and  Abaqus, 
independent  of  the  implementation  of  the  complex  perturbation  of  the  crack  length;  the  stress 
intensity  factors  calculated  by  CFEM  and  Abaqus  are  both  within  the  expected  errors  of  the 
reference  formulas  and  solution  methods  used.  As  described  in  Section  2.2.1,  the  imaginary  step 
size  has  an  effect  on  the  real-valued  solution  of  0(h2),  but  if  h  is  sufficiently  small,  the  J  integral 
calculation  will  not  be  significantly  effected  as  shown  in  this  analysis. 

Figure  17  shows  the  accuracy  in  K'CFEMIK'ref  as  a  function  of  the  J  integration  domain  for  the  4 

different  perturbation  methods.  Figure  18  shows  a  closer  view  of  the  same  data.  These  data 
indicate  from  which  domain  the  J'  integral  should  be  extracted  as  a  function  of  the  perturbation 
method.  From  Figure  17,  the  results  from  the  Tip  Node,  Domain  1,  and  Domain  4  methods  all 
converge  at  around  the  4th  integration  domain.  The  results  for  the  Face  Only  method  don’t 
converge  until  about  the  24th  integration  domain.  The  close  up  view  in  Figure  18  indicates  that 
ultimately,  the  best  results  are  obtained  from  the  36th  domain  using  the  Q4  method  to  implement 
the  imaginary  crack  length  perturbation. 
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Figure  17.  Imaginary  Displacement  Method  Influence  on  K'  Accuracy,  Infinite  Array, 
h/hx  =  10  6,  CFEM  K'  within  30%  of  Reference 


Figure  18.  Imaginary  Displacement  Method  Influence  on  K'  Accuracy,  Infinite  Array, 
h  !hx  =  10  6,  CFEM  K'  within  0.5%  of  Reference 


These  results  suggest  that  using  CTSE  to  estimate  the  first-order  sensitivity  of  K  with  respect  to  a 
change  in  a  is  more  accurate  when  larger  integration  domains  are  used  to  calculate  K.  Because  of 
this,  the  largest  7-integration  domain  (36)  is  used  to  examine  the  effects  of  varying  the  CTSE 
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complex  step  size  h  on  the  the  accuracy  of  K  and  K'  in  Figure  19.  The  results  show  the 
remarkable  stability  of  the  method  since  h!hx  varies  from  10  1  to  10  30°. 


Accuracy  of  K  and  K’  Varying  h 
(Q4  and  Tip  Only  Imaginary  Displacements) 

1  001000  1 - - - - - - 

1.000800  v - <; - <; - v. - <; - ;; - 

1.000600  - o  Kcfem/Kref  (ihQ4) - 

1.000400  - -  K'cfem/K'ref(ihQ4) 

1.000200  - a  Kcfem/Kref  (Tip  Only) - 

1.000000  - o  K'cfem/K'ref  (Tip  Only) - 

0.999800  I' - " -  |  oim.im.nm.in 
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Figure  19  Accuracy  of  K  and  IT  over  a  wide  range  of  Complex  Step  Size 
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3  Probabilistic  Sensitivity  Analysis  with  respect  to  Bounds  of 
Truncated  Distributions 


The  Score  Function  (SF)  method  is  a  mathematical  technique  to  compute  the  partial  derivative  of 
the  probability-of-failure,  mean  (juz)  or  standard  deviation  (cr2)  of  some  response  (Z)  with 
respect  to  the  parameters  of  the  input  PDFs,  e.g.  [34].  Under  this  program,  the  SF  method 

was  extended  to  address  the  case  of  computing  the  probabilistic  sensitivities  with  respect  to 
bounds  of  truncated  distributions.  Examples  of  truncated  distributions  include  uniform,  truncated 
normal,  and  truncated  Weibull  among  others. 

Bounds  on  variables  are  often  implemented  as  part  of  a  quality  control  program  to  ensure  a 
sufficient  pedigree  of  a  product  component  and  these  bounds  may  significantly  affect  the 
product’s  reliability  and  design  through  constraints  such  as  cost,  manufacturability  and 
reliability.  The  methodology  to  compute  the  sensitivity  of  the  probability-of-failure  with  respect 
to  bounds  of  truncated  distributions  is  summarized  below.  Extension  to  the  sensitivities  of  the 
mean  and  standard  deviation  of  the  response  can  be  found  in  ref.  [35] 

3.1  Methodology 

The  probability-of-failure  integral  can  be  written 


P,  =  J  ^(x)/x(x)</x 


(29) 


where  x  represents  a  vector  of  random  variables,  /(x)  denotes  the  indicator  function,  (1  if  failure 
occurs,  0  otherwise),  and  fx  represents  the  joint  probability  density  function  of  the  random 
variables. 

The  derivative  of  the  probability  integral  with  respect  to  a  parameter  of  a  random  variable  that 
affects  the  boundary  can  be  determined  by  using  the  idea  of  the  classical  Material  Derivative 
from  continuum  mechanics  [36], 

£  |  W(x,t)dV  =  |  dV+  |  x¥(x,t)vJnJdS  (30) 

where  T  is  a  property  of  the  continuum,  y  denotes  the  velocity  of  the  material,  n  represents  the 
unit  normal  along  the  boundary  S,  and  V  is  the  volume  enclosed  by  S.  The  surface  integral  term 
in  Equation  30  is  the  value  of  T  on  the  boundary  multiplied  by  the  volume  swept  by  the  particles 
on  the  boundary  in  the  time  interval  dt,  integrated  over  dS.  This  term  can  be  considered  as  a 
flux  of  the  property  lT  over  the  surface  S. 
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The  concept  of  the  material  derivative  can  be  utilized  to  take  the  derivative  of  the  probability 
integral,  Equation  30,  with  respect  to  a  bound  of  a  random  variable  PDF.  Here,  the  independent 
parameter,  0r  is  a  bound  of  the  distribution  representing  X,  the  JPDF  is  equivalent  to  the 
property  T,  the  volume  is  the  N  dimensional  random  variable  space,  and  S  is  the  surface  of  the 
random  variable  space  remaining  when  random  variable  X.  is  set  to  the  bound  6r  For 
rectangular  truncation,  the  surface  S  is  straightforward  to  compute  as  the  independent  parameter, 
0i,  is  a  bound  of  the  N  dimensional  random  variable  space. 

The  unit  normal  and  the  equivalent  velocity  term  and  their  relation  can  be  discerned  from  a 
problem  of  two  random  variables,  see  Figure  20.  Since  the  independent  parameter  0  is  an 
element  of  X,  the  velocity  becomes  v  =  ckl '36  =  1.  At  the  lower  bound  v  and  n  are  in  opposite 
directions,  hence,  the  dot  product  v  ./z .  equals  -1.  At  the  upper  bound  vjn  j  equals  +1. 


Figure  20  Descriptions  of  Velocity  and  Unit  Normal  Along  Bounds 

Carrying  out  the  Material  Derivative  of  the  probability-of-failure  integral,  the  sensitivities 
become 


3Pf 


=  fXi(a)(Pf-Pn 


(31) 


3Pf 

cbt 


=-fXimpf-p f) 


(32) 


where  P,  represents  the  usual  probability-of-failure  integral  and  P“  and  P/  denote  supplemental 
flux  integrals.  Thus,  the  effort  to  obtain  the  sensitivities  becomes  one  of  evaluating  the  flux 
integrals  P“  and  P/  in  addition  to  the  probability-of-failure. 
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The  flux  integral  can  be  evaluated  by  integrating  the  joint  probability  density  function  times  the 
indicator  function  over  the  surface  S  defined  by  the  condition  x,  =  6.  That  is,  one  integrates  the 
JPDF  in  the  failure  region  over  the  surface  S.  Thus,  the  dimension  of  the  flux  integral  is  N-l.  A 
schematic  of  the  flux  in  two  dimensions  (two  random  variable  problem)  is  shown  in  Figure  21 
where  the  value  of  the  joint  PDF  along  the  upper  bound  of  Xx  is  outlined. 


Figure  21  Flux  of  the  JPDF  in  the  Failure  Region  over  the  Surface  S 

It  is  straightforward  to  calculate  the  flux  integrals  using  sampling  as 


1 

m 


k=l 


(33) 


That  is,  the  samples  used  to  compute  Pf  are  reused  to  compute  Pf  by  projecting  the  samples 
onto  the  surface  Xt  =  a  or  A  =  b  and  the  indicator  function  is  evaluated  along  the  surface 

x.  =  er 


3.2  Results  and  Discussion 

3.2.1  Numerical  Example 

A  two  dimensional  example  is  solved  to  illustrate  the  concepts.  Consider  a  limit  state  of 
g(r,s )  =  r—s  with  R  and  S  modeled  as  independent  random  variables.  The  indicator  function 
defines  the  failure  region  as 
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l[r,  ,v]  =  1  if  g(r,  .v)  <  0 

0  otherwise 


(34) 


R  is  a  standard  normal  distribution  and  S  is  a  uniform  distribution  with  bounds  a  =  0  and  b  =  1. 


1 

fR  (r)  =  , —  Exp\-r  / 2]  -  oo  <  r  <  oo 

•v2;r 

/s(*)=l  0<s<l 


0 


otherwise 


(35) 

(36) 


Table  17  summarizes  results  obtained  using  both  the  flux-based  methodology  and  the  standard 
finite  difference  (forward  differencing)  method.  The  derivatives  of  the  probability-of-failure  with 
respect  to  the  bounds  were  computed  using  Eqs.  (31)  and  (32)  using  both  symbolic  integration 
and  sampling.  Finite  difference  estimates  were  also  computing  using  symbolic  integration  and 
sampling  with  a  forward  step  size  of  0.00001. 


Table  17  Probabilistic  Sensitivity  Results  for  Limit  State  g  =  r-s  (R  Standard  Normal,  S 

Uniform) 


Flux-based 

Finite  Difference 

Integration 

Sampling 

Integration 

Sampling 

(exact) 

(104  samples) 

(exact) 

(106  samples) 

dPf 

0.1844 

0.18421 

0.1844 

.1876' 

da 

dPf 

0.1560 

0. 157 1 1 

0.1570 

On 

O 

OO 

cb 

'expected  value  of  100  trials 


The  sensitivities  for  the  flux-based  and  finite  difference  methods  using  exact  integration  are  very 
close,  as  expected.  The  results  using  sampling  are  also  in  good  agreement;  however,  note  that  104 
samples  were  used  for  the  flux-based  approach  to  obtain  a  solution  with  good  accuracy  versus 
106  samples  for  the  finite  difference  sampling-based  approach.  The  superiority  of  the  flux-based 
approach  compared  to  the  standard  finite  difference  (forward  differencing)  approach  using 
sampling  can  be  shown  clearly  by  comparing  the  95%  confidence  bounds  and  the  coefficient  of 
variation  (COV  =  standard  deviation/mean)  obtained  using  both  methods  obtained  from  100 
trials.  Figure  22  shows  a  plot  of  the  95%  confidence  bounds  for  finite  difference  (dashed)  versus 
flux-based  (solid)  for  dP,  lea  (similar  results  are  obtained  for  cPf  Icb).  The  number  of  samples 
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are  shown  in  the  X  axis  in  log10  scale.  The  bounds  for  the  flux-based  approach  are  so  much 
narrower  than  the  finite  difference  approach  that  they  show  almost  as  a  straight  line.  The  bounds 
are  so  wide  for  the  finite  difference  method  that  any  solution  obtained  is  completely  unreliable 
until  the  number  of  samples  approaches  one  million. 

A  closer  examination  of  the  benefits  of  using  common  random  variables  during  the  sampling 
process  is  provided  in  Table  18.  Variance  results  for  negative,  approximately  zero,  and  positive 
sampling  correlation  are  provided.  The  results  clearly  show  that  positive  correlation  provides 
approximately  a  three  times  reduction  in  the  variance  of  the  sensitivities  with  respect  to 
independent  sampling;  at  zero  cost.  Positive  correlation  was  accomplished  simply  by  using  the 
same  samples  to  compute  the  probability-of-failure  and  the  fluxes. 
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Table  18  Variance  of  sensitivity  estimates  dPf  IdO  as  a  function  of  sampling  correlation 
(100  trials)  -  actual  correlation  in  parentheses 


# 

dPf 

dPf 

dPf 

dPf 

dPf 

dPf 

Samples 

da 

da 

da 

db 

cb 

db 

(p<  0) 

(P-0) 

(p>  0) 

(P<  0) 

(P-0) 

(p>  0) 

103 

8.81E-4 

(p=-0.75) 

4.64E-4 

(P~0) 

1.37E-4 

(p=0.82) 

4.84E-4 

3.48E-4 

(P~0) 

1.29E-4 

(p=0.78) 

(p=-0.46) 

104 

7.37E-4 

4.01E-5 

1.44E-5 

4.03E-4 

2.70E-5 

1.22E-5 

(p=-0.66) 

(p=0.09) 

(p=0.68) 

(p=-0.27) 

(p=0.09) 

(p=0.64) 

Table  19  shows  the  coefficient  of  variation  of  the  two  methods  as  a  function  of  the  number  of 
samples.  The  COV  for  the  finite  difference  method  is  approximately  two  orders  of  magnitude 
larger  than  that  obtained  using  the  flux-based  method  for  the  same  number  of  samples.  These 
results  imply  that  there  are  approximately  4  orders  of  magnitude  difference  in  the  number  of 
samples  required  to  achieve  similar  accuracy  using  both  methods.  The  explanation  is  that  the 
finite  difference  method  requires  an  approximation  of  a  limiting  process  estimated  by  subtracting 
two  near-equal  numbers.  The  flux-based  approach,  on  the  other  hand,  requires  no  limiting 
process  or  subtraction  of  near-equal  numbers.  The  results  from  Table  19  indicate  that  using  the 
flux-based  approach  with  103  samples  is  superior  to  the  finite  difference  method  with  106 
samples. 
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Table  19  Coefficient  of  Variation  (100  Trials)  of  Sensitivities  with  respect  to  Bounds  as  a 

Function  of  the  Number  of  Samples 


Flux-based 

Finite  Difference 

#  Samples 

cPf 

dPf 

dPf 

dPf 

da 

cb 

da 

cb 

103 

.076 

.063 

4.4 

8.9 

104 

.019 

.022 

2.1 

1.8 

105 

.0076 

.0076 

.61 

0.66 

106 

.0021 

.0025 

.16 

0.22 
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4.  Conclusions 


Development  of  probabilistic  sensitivities  is  frequently  considered  an  essential  component  of  a 
probabilistic  analysis  and  often  critical  towards  understanding  the  physical  mechanisms 
underlying  failure  and  modifying  the  design  to  mitigate  and  manage  risk.  As  a  result,  it  is  useful 
to  explore  and  enhance  the  probabilistic  sensitivity  methods  that  are  available  for  analysis.  In  this 
work,  three  sensitivity  approaches  were  developed  and  applied.  In  particular,  the  developments 
accomplished  include:  1)  a  new  method  to  compute  the  sensitivity  of  the  probability-of-failure 
with  respect  to  the  parameters  that  define  the  probability  of  detection  curve  of  an  NDE,  2)  the 
investigation  of  complex  variable  methods  for  sensitivity  analysis  of  structures  using  finite 
element  analysis,  and  3)  enhancements  to  the  score  function  method  to  compute  the  sensitivity 
with  respect  to  bounds  of  truncated  distributions. 


4.1  Conclusions 

A  probabilistic  methodology  was  developed  and  verified  to  compute  the  sensitivity  of  the 
probability-of-failure  with  respect  to  parameters  of  a  POD  curve  used  during  an  inspection 
process.  The  methodology  has  several  attractive  features.  The  formulation  is  such  that  the 
probabilistic  sensitivities  can  be  obtained  at  negligible  cost  using  standard  (Monte  Carlo)  or 
other  sampling  (Latin  hypercube,  importance  sampling)  methods.  The  key  is  that  the  same 
samples  used  to  compute  probability-of-failure  are  reused  to  compute  the  sensitivities.  As  such, 
the  methodology  can  be  implemented  in  a  post-processing  non-intrusive  manner  thereby 
facilitating  implementation  with  existing  or  commercial  codes;  the  methodology  only  needs  the 
crack  sizes  at  the  times  of  inspection  and  the  cycles-to-failure  for  each  realization. 

The  formulation  is  generic  and  not  limited  to  any  specific  random  variables,  fracture  mechanics 
formulation,  or  any  specific  POD  curve  as  long  as  the  POD  is  modeled  parametrically.  The 
example  problems  demonstrated  sensitivities  with  respect  to  two  parameters  of  the  lognormal- 
form  POD  curve  but  any  POD  form  with  any  number  of  parameters  can  be  used  if  the  Q  term 
can  be  calculated. 

The  resulting  equations  are  remarkably  simple  and  involve  an  evaluation  and  summation  of  an  Q 
term  with  the  summation  occurring  only  for  samples  that  have  failed  before  time  t  and  the 
evaluation  occurs  only  at  the  time  of  inspection.  The  form  of  Q  requires  the  derivative  of  the 
POD  curve  with  respect  to  its  parameters  and  can  be  derived  analytically.  The  time  dependent 
nature  of  the  sensitivity  is  accounted  for  by  the  time  dependent  indicator  function  term.  The 
sensitivity  equations  for  a  particular  inspection  are  unchanged  in  form  due  to  subsequent 
inspections  however  their  numerical  evaluation  may  change  due  to  samples  being  detected  and 
removed  by  subsequent  inspections. 

4.2  Conclusions 

Complex  Taylor  series  expansion  (CTSE),  Fourier  differentiation  (FD),  and  central  differencing 
(CD)  can  all  be  used  to  calculate  shape  sensitivities.  To  our  knowledge,  the  work  presented  here 
marks  the  first  time  that  FD  has  been  used  for  calculating  finite  element  shape  sensitivities. 
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Experience  shows  that  each  method  can  compute  reasonably  accurate  first  order  derivatives; 
however,  the  complex  variable  methods  are  attractive  in  that  they  are  extremely  robust  to  the 
complex  step  size  h  whereas  finite  difference  is  not.  For  2-D  finite  element  problems  using 
second  order  shape  functions,  the  error  in  the  model  was  greater  than  the  error  due  to  the 
truncation  errors  associated  with  the  numerical  differentiation  methods.  This  means  that  the 
complex  variable  sensitivity  methods  do  not  offer  extra  accuracy  compared  to  CD  if  the  model  is 
not  sufficiently  accurate.  They  do,  however,  offer  some  ease  in  implementation.  If  the  model 
were  made  to  be  more  accurate,  such  as  with  higher  order  basis  functions  or  more  elements,  then 
FD  and  CTSE  could  offer  improved  accuracy  with  respect  to  CD. 

It  has  been  shown  that  FD  is  capable  of  producing  highly  accurate  sensitivities  for  functions  with 
solutions  that  are  accurate  to  machine  precision  [27].  The  trade-off  for  the  increased  accuracy  of 
FD  is  the  requirement  of  several  complex  sample  points.  It  may  take  as  much  as  three  times  more 
computational  effort  to  generate  a  complex  sample  than  a  real  valued  sample.  Thus  for  the 
problems  described  in  this  paper,  FD  is  not  a  good  choice  for  shape  sensitivity  calculations,  due 
to  the  limited  accuracy  of  the  solution.  CTSE  requires  only  half  of  the  number  of  sample  points 
required  by  CD,  thus  CTSE  only  requires  about  1.5  times  more  computational  effort  than  CD,  at 
most.  When  the  evaluation  of  a  shape  sensitivity  requires  displacement  of  only  a  small 
proportion  of  nodal  ordinates  the  additional  computational  effort  may  be  much  less.  This  coupled 
with  the  fact  that  CTSE  doesn’t  require  changing  the  location  of  nodes  in  the  complex  plane 
means  that  it  may  still  be  a  good  choice  for  shape  sensitivity  problems. 

One  of  the  biggest  problems  with  CD  is  that  it  requires  the  user  to  change  the  location  of 
multiple  nodes  and  maybe  surrounding  elements  to  implement  a  change  in  the  desired  shape. 

This  may  lead  to  elements  with  poor  aspect  ratios,  or  require  complete  remeshing  of  the  domain, 
or  both.  This  is  especially  true  when  the  step  size  is  rather  large.  Since  CTSE  does  not  require 
moving  the  nodes  in  the  real  plane,  remeshing  is  not  required,  and  there  is  less  concern  over  the 
aspect  ratio  of  the  elements.  Furthermore,  if  only  the  first  order  sensitivity  is  required,  than  the 
step  size  can  be  made  very  small  without  fear  of  increasing  the  round-off  error. 

4.3  Conclusions 

Efficient  evaluation  of  the  sensitivity  of  the  probability-of-failure  or  the  response  moments  to  the 
bounds  of  truncated  distributions  can  provide  useful  information  in  the  design  stage  in  order  to 
optimize  product  reliability,  minimize  cost,  determine  quality  assurance  procedures,  etc.  The 
method  outlined  here  can  be  used  to  compute  these  sensitivities  with  a  significant  improvement 
in  computational  efficiency  over  standard  finite  difference  methods. 

The  methodology  is  based  upon  an  application  of  the  material  derivative  concept  to  the 
probability-of-failure  or  the  response  moment  integrals  thus  yielding  a  flux  integral  that  must  be 
computed  in  addition  to  the  standard  probability-of-failure  or  response  integrals.  The 
methodology  is  applicable  to  any  limit  state  formulation,  either  component  or  system,  and  any 
random  variables  described  by  a  truncated  joint  probability  density  function  containing  either 
correlated,  e.g.,  truncated  multivariate  normal,  or  independent  random  variables. 

The  sensitivities  require  a  supplemental  flux  integral  for  each  bound  that,  when  combined  with 
the  probability-of-failure  and  kernel  functions,  provides  the  sensitivity  with  respect  to  the  bound 
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of  a  truncated  distribution.  However,  the  flux  integral  itself  is  a  probability  integral  and, 
therefore,  amenable  to  solution  using  existing  probabilistic  methods.  If  sampling  is  used  to 
compute  the  probability-of-failure,  the  samples  can  be  reused  to  compute  flux  integrals  by 
projecting  the  samples  to  the  bound  of  interest.  Thus,  the  flux  integrals  are  computed  with 
negligible  additional  computation. 

The  superiority  of  the  flux-based  approach  over  the  standard  finite  difference  method  was  clearly 
evident  from  numerical  studies  using  Monte  Carlo  sampling  that  indicated  that  the  estimate  of 
the  sensitivities  using  the  flux-integral  approach  required  approximately  4  orders  of  magnitude 
fewer  samples  for  the  same  accuracy  as  a  standard  finite  difference  approach. 
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5.  Recommendations 

The  field  of  probabilistic  sensitivity  analysis  is  a  critical  aspect  of  probabilistic  analysis.  In  order 
to  put  the  results  of  a  probabilistic  analysis  into  practice,  it  is  necessary  to  know  the  relative 
importance  of  the  variables  to  the  problem  at  hand.  This  is  very  dynamic  and  evolving  field  with 
the  work  describing  only  portion  of  the  methods  in  the  literature.  In  effect,  the  field  largely  did 
not  exist  20  years  ago.  Hence  it  is  incumbent  upon  the  analyst  to  stay  abreast  of  the  latest 
developments. 

Although  significant  progress  has  been  made,  there  are  still  challenges  to  be  overcome.  The 
tension  of  computational  efficiency  versus  accuracy  is  a  continual  challenge.  For  example,  the 
Score  Function  method  provides  sensitivities  for  negligible  cost  since  the  Monte  Carlo  samples 
are  reused,  however,  the  variance  of  the  sensitivities  may  be  too  large.  As  a  result,  accurate 
quantification  of  the  sensitivities  may  require  significantly  more  samples  than  used  to  compute 
the  probability-of-failure  or  the  response  moments.  This  is  likely  to  be  true  for  variables  of 
secondary  importance  but  maybe  true  for  primary  variables  also.  Thus,  new  methods  that  provide 
accurate  sensitivity  information  for  minimal  cost  are  still  needed. 

The  complex  variable  sensitivity  methods  investigated  provide  an  exciting  new  approach  to 
calculate  sensitivities  of  all  types.  This  is  especially  true  for  shape  sensitivities  as  the  only  the 
affected  nodes  need  to  be  perturbed  along  the  imaginary  axis  -  no  modifications  to  the  real  mesh 
are  required.  Since  the  mathematics  are  general,  it  is  incumbent  that  we  explore  further 
applications  and  derivations  of  this  method. 
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LIST  OF  ACRONYMS,  ABBREVIATIONS,  AND  SYMBOLS 

ACRONYM  DESCRIPTION 

POD  Probability-of-Detection 

Pf  Probability-of-failure 

fx  Joint  density  function 

I  Indicator  function 

X  Vector  of  random  variables 

a  Crack  size 

t  Time  or  cycles-to-failure 

0  Parameter  of  POD  curve 

Q  Kernel  function  for  POD  sensitivity  analysis 

C  Paris  crack  growth  constant 

/ n  Paris  crack  growth  exponent 

ac  Critical  crack  size 

g  Limit  state 

K1C  Fracture  toughness 

K  Stress  intensity  factor 

K'  Derivative  of  K  with  respect  to  crack  length 
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